Esercitazione 4

Compensazione dinamica

Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.

1 Funzioni di utilità

# Generazione sinusoidale con diverse armoniche
signal <- function(t, pars, rad = FALSE) { 
  stopifnot(is.data.frame(pars))
  with(pars, {
    if (!rad) {
      phi <- phi/180*pi
      f <- 2*pi*f
    }
    map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
  })
}

2 Simulazione uscita sistema di misura

# Parametri del misurando (input del sistema di misura)
pars <- tibble(
  w = c(1, 0.3, 0.3), 
  f = c(1/30, 1/10, 1/4), 
  phi = c(0, 0, 0)
)

# Generiamo il misurando u
Ts <- 0.01            # periodo di campionamento
fs <- 1/Ts            # freq. di campionamento   
u <- tibble(
  t = 0:10000 * Ts,
  u = signal(t, pars),
  un = u + rnorm(length(t), 0, pars$w[1]/10)
)

# Definiamo la Funzione di Trasferimento dello strumento
num <- c(5, 1)   # 5s + 1
den <- c(1, 2, 1) # s^2 + 2s + 1
H <- tf(num, den)
print(H)

y1:
    5 s^1 + 1 
  - - - - - - - - - -
     s^2 + 2 s + 1 


Transfer Function: Continuous time model 

# Simulazione dell'uscita nel tempo
tmp <- u %>% mutate(
    out = lsim(H, un, t)$y[1,],
    ideal = lsim(H, u, t)$y[1,]
  ) %>% mutate(n = row_number(), .before = t)

tmp %>% plot_ly(x = ~t) %>%
  add_lines(y = ~un, name = "Input") %>%
  add_lines(y = ~ideal, name = "Output ideale") %>%
  add_lines(y = ~out, name = "Output") %>%
  layout(xaxis=list(title="Tempo (s)"), yaxis=list(title=""), title = "Simulazione dell'uscita")

Notare che lsim calcola la risposta nel tempo di un sistema lineare.

2.1 Spettro dei segnali misurati

Ta <- max(u$t)
un.fft <- compute_fft(tmp, Ta, un)
ideal.fft <- compute_fft(tmp, Ta, ideal)
out.fft <- compute_fft(tmp, Ta, out)

bind_rows(
  un.fft %>% mutate(caso = "Input"),
  ideal.fft %>% mutate(caso = "Output ideale"),
  out.fft %>% mutate(caso = "Output")
) %>% dplyr::filter(f < 1) %>% 
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#F8766D") +
  geom_point(colour = "#F8766D") +
  facet_wrap(~ caso, ncol = 1, scales = "free_x") +
  labs(
    x = "Frequenza (Hz)",
    y = "Intensità"
  )

NB: non si vede, ma il rumore bianco è presente. Provare ad ingrandire la dev. std. del rumore addizionato a un.

3 Stima del misurando a partire dal segnale di uscita

# Segnale in uscita con rumore
yI <- lsim(H, u$u, u$t)$y[1,]
yn <- yI + rnorm(length(u$t), 0, pars$w[1]/10)


# Ricaviamo l'inversa della Funzione di Trasferimento:
Hinv <- tf(den, num)

[1] "TFCHK: Transfer function may not be proper and may lead to errors. Num > Den"

print(Hinv)

y1:
     s^2 + 2 s + 1 
  - - - - - - - - - -
    5 s^1 + 1 

Transfer Function: Continuous time model 

# Simulazione dell'ingresso
#input <- lsim(Hinv, yn, y$t, x0 = rep(0, length(pole(Hinv))))

Error in tf2ss(num, den) :
TF2SS: The order of the Numerator should be equal or lesser than the Denominator.

Implementazione della compensazione dinamica.

Il problema risiede nel fatto che la funzione lsim è progettata per calcolare la risposta nel tempo di un sistema LTI proprio (\(deg(\text{NUM}) \leq deg(\text{DEN})\)). Nel caso dell’inversa \(H^{-1}(j\omega)\), è come se si stimasse la risposta di un sistema improprio, per cui lsim fallisce.

Risolvo elaborando tutto nel dom della freq e poi antitrasformo.

# filtro l'uscita per ridurre il rumore
y <- u %>% select(t) %>% mutate(n = row_number(t), .before = t) %>%
  mutate(
    y_id = yI,
    y_n = yn
  )

compute_fft(y, Ta, yn) %>% dplyr::filter(f < 2.5) %>%
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#56B4E9") +
  geom_point(colour = "#56B4E9") +
  labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro dell'uscita con rumore")

# filtro a 1Hz
fc <- 1
flt <- butter(2, fc/(fs/2))

ggbodeplot_digital(flt, fs, fmin = 1e-1, fmax = 1e1) + geom_vline(xintercept = fc, color="red", linetype = 2)

y <- y %>% mutate(y_flt = signal::filtfilt(flt, yn))

# per padding automatico di filtfilt ho una distorsione, perciò taglio il segnale di uscita
y <- y %>% dplyr::filter(t <= 99)

# Nuovo tempo di campionamento
Ta <- max(y$t)

# Windowing perché segnale complessivo non periodico -> FFT non funziona
N <- length(y$t)
y <- y %>% mutate(
  win = signal::hanning(N),
  y_win = y_flt * win
)

y %>% pivot_longer(cols = c("y_win", "y_flt")) %>% 
  ggplot(aes(x = t, y = value, colour = name)) +
  geom_line() +
  geom_line(aes(y=win), color = "red", linetype = 22) +
  geom_line(aes(y=-win), color = "red", linetype = 22) +
  scale_color_discrete(labels = c(y_win = "Windowing", y_flt = "Filtrata")) +
  labs(x = "Tempo (s)", y="Output", colour = "Segnale")

[ChatGPT] Per chiarezza e per evitare errori numerici, viene calcolata esplicitamente la funzione di risposta in frequenza, anziche usare la freqresp.

# Trasformata di Fourier (FFT)
y <- y %>% mutate(
  Y =fft(y_flt),
  Y_win = fft(y_win)
)

# Asse frequenze
freq <- tibble(
  k = 0:(N-1),
  f = ifelse(k <= N/2, k, k-N) * fs/N,   # FFT simmetrica nell'asse delle freq.
  w = 2*pi*f
)

# Calcolo della funz. di risposta in frequenza
Hiw <- freqresp(H, freq$w)

# Inverto la caratteristica dinamica e anti-trasformo
u_comp <- tibble(
  t = y$t,
  U = y$Y / Hiw,
  U_win = y$Y_win / Hiw,
  u = Re(fft(U, inverse = TRUE)) / N,
  u_win = Re(fft(U_win, inverse = TRUE)) / N # Re() serve a trascurare dei residui di calcolo nella parte Im 
  # la fft(inverse = TRUE) ritorna l'inversa NON-normalizzata, perciò si divide per la lunghezza del campione
)

u_comp %>% plot_ly(x = ~t) %>% #add_lines(y = y$y_n, name = "Uscita") %>%
  add_lines(y = u$un[u$t <= Ta], name = "Reale") %>%
  add_lines(y = u$u[u$t <= Ta], name = "Ideale") %>%
  add_lines(y = ~u, name = "Compensato") %>%
  add_lines(y = ~u_win, name = "Windowing")

3.1 Provo con padding simmetri a zero [ChatGPT]

# Ampiezza padding
M <- 4*N

pad_left  <- floor((M - N)/2)
pad_right <- ceiling((M - N)/2)

t_start <- y$t[1] - pad_left * Ts
y_ext <- tibble(
  t = seq(from = t_start, by = Ts, length.out = M),
  y_pad = c(
    rep(0, pad_left),
    y$y_flt,
    rep(0, pad_right)
  )
)

# Verifica che il segnale sia centrato: deve restituire TRUE
all.equal(
  y_ext$y_pad[(pad_left+1):(pad_left+N)],
  y$y_flt
)

[1] TRUE

y_ext %>%
  ggplot(aes(t, y_pad)) +
  geom_line() +
  geom_line(data = y, aes(x = t, y = y_n), color = "red")

freq <- tibble(
  k = 0:(M-1),
  f = ifelse(k <= M/2, k, k-M) * fs/M,   # FFT simmetrica nell'asse delle freq.
  w = 2*pi*f
)

Hiw <- freqresp(H, freq$w)

y_ext <- y_ext %>% mutate(
  Y_pad = fft(y_pad),
  U_pad = Y_pad/Hiw,
  u_pad = Re(fft(U_pad, inverse = TRUE)) / M
)

u_comp <- u_comp %>% mutate(
  u_pad = y_ext$u_pad[(pad_left+1):(pad_left+N)]
)

u_comp %>% plot_ly(x = ~t) %>% #add_lines(y = y$y_n, name = "Uscita") %>%
  add_lines(y = u$un[u$t <= Ta], name = "Reale") %>%
  add_lines(y = u$u[u$t <= Ta], name = "Ideale") %>%
  add_lines(y = ~u, name = "Compensato") %>%
  add_lines(y = ~u_win, name = "Windowing") %>%
  add_lines(y = ~u_pad, name = "Padded")

Applico un filtro sul segnale padded:

compute_fft(u_comp, Ta, u_pad) %>% dplyr::filter(f < 2) %>%
  ggplot(aes(x = f, y = mod)) +
  geom_spoke(aes(angle = pi/2, y=0, radius = mod), linewidth = 0.2, colour = "#56B4E9") +
  geom_point(colour = "#56B4E9") +
  labs(x = "Frequenza (Hz)", y = "Intensità", title = "Spettro dell'ingresso compensato")

fc <- 2
flt <- butter(2, 2*fc/fs)

u_comp <- u_comp %>% mutate(
  u_flt = signal::filtfilt(flt, u_pad)
)

u_comp %>% plot_ly(x = ~t) %>% 
  add_lines(y = u$un[u$t <= Ta], name = "Reale") %>%
  add_lines(y = u$u[u$t <= Ta], name = "Ideale") %>%
  add_lines(y = ~u_pad, name = "Padded") %>%
  add_lines(y = ~u_flt, name = "Filtered")

ora accetto il risultato